load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"
load "/public/users/wrf/Codes/WRF_scripts/utils/substring.ncl"

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;            User Config                   ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

if (.not. isvar("configFile")) then
    configFile = "/public/users/wrf/Codes/WRF_scripts/WRF_Post/StationLists/sswdn_stnlist.csv"
    print("configFile is undefined, using default: " + configFile)
else
    print("Using configFile: " + configFile)
end if

lines = readAsciiTable(configFile,1,"string",1)
delim = ","

stnnum    = ndtooned(str_get_field(lines,1,delim))
stnlist    = ndtooned(str_get_field(lines,2,delim))
stnlon    = ndtooned(tofloat(str_get_field(lines,3,delim)))
stnlat    = ndtooned(tofloat(str_get_field(lines,4,delim)))
stndomain = ndtooned(str_get_field(lines,5,delim))
nstns     = dimsizes(stnlist)

; default output time interval is 30min
if (.not. isvar("output_interval")) then
    output_interval = 0.25 ; 15 min
    print("output_interval is undefined, using default = " + output_interval)
else
    print("Using output_interval = " + output_interval)
end if

; test if input_file is given
if (.not. isvar("input_file")) then
    input_file = "/public/users/wrf/wrfout/2012103120/wrfout_d02_2012-10-31_12:00:00"
    print("input_file is undefined, using default = " + input_file)
else
    print("Using input_file: : " + input_file)
end if

filename = systemfunc("basename " + input_file)

; test if output_prefix is given
if (.not. isvar("output_prefix ")) then
    output_prefix = "~/wrfpost/test/station_sswdn"
    print("output_prefix is undefined, using default = " + output_prefix)
else
    print("Using output_prefix: " + output_prefix)
end if

if (.not. isfilepresent(output_prefix)) then
    system("mkdir " + output_prefix)
end if

; check start_date
if (.not. isvar("start_date")) then
    start_date = substring(filename, 11, 23)
    print("start_date is undefined, using default = " + start_date)
else
    print("Using start_date: " + start_date)
end if

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;           Open NetCDF File               ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
a = addfile(input_file + ".nc", "r")

domain = substring(filename, 7, 9)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;              Process Time                ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
times = a->Times
hours = WRF_Times2Udunits_c(times, 0)
wrf_dt = hours(1) - hours(0)
; UTC -> CST
hours = hours + 8 ; 
timeStr_cst = ut_string(hours, "%Y-%N-%D_%H:%M:%S")
ntimes = dimsizes(hours)

dt = tointeger(output_interval / wrf_dt)

; create time series in desired time resolution
output_hours = fspan(hours(0), hours(ntimes-1), tointeger((hours(ntimes-1) - hours(0))/output_interval + 1))
output_ntime = dimsizes(output_hours)
copy_VarAtts(hours, output_hours)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;               Read SSWDN                 ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
sswdn = wrf_user_getvar(a,"SSWDN",-1)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;          Loop for each station           ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

do istn = 0, nstns - 1, 1
    ; if not in the desired domain
    if (.not. isStrSubset(filename, stndomain(istn))) then
        continue
    end if

    ; Setting Output Filename
    output_filename = output_prefix + "/sswdn_output_" + stnnum(istn) + "_" + stnlist(istn) +  "_" + domain +  "_" + ut_string(hours(0), "%Y-%N-%D_%H") + ".txt"

    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ;             Lon/Lat to XY                ;
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    res = True
    res@returnInt = False   ; return real values
    loc = wrf_user_ll_to_ij(a, stnlon(istn), stnlat(istn), res) - 1

    sxc0 = floattointeger(loc(0))
    syc0 = floattointeger(loc(1))

    syc1 = syc0 + 1
    sxc1 = sxc0 + 1

    yy = (/syc0, syc1/)
    xx = (/sxc0, sxc1/)
    kk = (/0, 1, 2/)

    print ("X/Y location is: " + loc)

    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ;                   read SSWDN                     ;
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    tmp_sswdn = rm_single_dims(linint2( xx, yy, wrf_dyy_getvar(a,"SSWDN",-1, kk, yy, xx), False, loc(0), loc(1), 0))
    sswdn = dyy_interp_ts(hours, tmp_sswdn, output_hours, "linear")
    delete(tmp_sswdn)

    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ; output precipitation tendency for dt hours       ;
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ; station number and time
    header = "ID,  SITE NAME,     TIME"
    printArray = stnnum(istn) + "," + stnlist(istn) + "," + ut_string(output_hours, "%Y-%N-%D_%H:%M:%S")

    ; SSWDN
    header = header + ", sswdn" 
    printArray = printArray + ", " + sprintf("%7.2f",sswdn)

    ; combine header and body
    allPrintArray = array_append_record(header, printArray, 0)
    asciiwrite(output_filename, allPrintArray)
    print("Creating: " + output_filename)

    ; delete temporary variables
    delete(header)
    delete(printArray)
    delete(allPrintArray)
end do ; loop for station
